The emergence of coherence in complex networks of heterogeneous dynamical systems 
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We present a general theory for the onset of coherence in collections of heterogeneous maps 
interacting via a complex connection network. Our method allows the dynamics of the individual 
uncoupled systems to be either chaotic or periodic, and applies generally to networks for which the 
number of connections per node is large. We find that the critical coupling strength at which a 
transition to synchrony takes place depends separately on the dynamics of the individual uncoupled 
systems and on the largest eigenvalue of the adjacency matrix of the coupling network. Our theory 
directly generalizes the Kuramoto model of equal strength, all-to-all coupled phase oscillators to the 
case of oscillators with more realistic dynamics coupled via a large heterogeneous network. 

PACS numbers: 05.45.-a, 05.45.Xt, 89.75.-k 



In recent years, much progress has been made in de- 
scribing the complex structure of real world networks 
0, @- The study of dynamical processes taking place 
in such complex networks has applications in fields rang- 
ing from biology to engineering. One of the most impor- 
tant phenomena involving networks of coupled dynamical 
systems is the emergence of large-scale coherent behavior 
H,Q- It is often observed that large collections of hetero- 
geneous dynamical systems (e.g., cells, fireflies) synchro- 
nize their rhythms so that a significant proportion of the 
systems have states that are highly correlated with those 
of the others. It is natural to ask what determines the 
emergence of such coordinated behavior given the net- 
work of interactions between the dynamical systems and 
their individual dynamics. 

The case of equal-strength, all-to-all coupled phase os- 
cillators was studied by Kuramoto who considered 
the case of N oscillators, each of which is described by 
a phase 9j and an intrinsic frequency ojj . Kuramoto as- 
sumed sinusoidal coupling so that the phase of oscillator 
j evolves as 6j = cjj + fc52 m =i sm (^m — Kuramoto 
found that, in the limit N —> oo, for coupling strengths k 
less than a critical coupling strength k c that depends on 
the distribution of frequencies, the phases of the oscilla- 
tors are incoherent, i.e., Oj are uniformly distributed on 
[0, 2ir). For values of the coupling strength k larger than 
fe c , a significant fraction of the oscillators evolve with a 
common frequency. The Kuramoto model has become a 
classical paradigm for the emergence of coherent behav- 
ior in an ensemble of heterogeneous oscillators (see jg] for 
reviews) . 

Although the Kuramoto model has the advantage that 
it can be treated analytically, it depends on the sys- 
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tern being simple in two major aspects. First, the net- 
work is assumed to be all-to-all, so that every oscillator 
is coupled with uniform strength to every other oscilla- 
tor. Second, the dynamics and coupling term are highly 
simplified: each oscillator is described only by its phase 
and the coupling term is sinusoidal. Recently, collec- 
tions of coupled dynamical systems which have either 
a more general interaction network or a more general 
dynamics have been studied. In previous works 0, Q 
(see also H E3, [H E E3), we have studied the Ku- 
ramoto phase oscillator model generalized to the case of 
a general interaction network described by an adjacency 
matrix. We found that there is still a transition to syn- 
chrony at a critical coupling strength that depends on 
the largest eigenvalue of the adjacency matrix and the 
distribution of frequencies. On the other hand, 'globally' 
coupled (i.e., equal coupling strength, all-to-all) collec- 
tions of many dynamical systems with more general dy- 
namics (e.g., mixed collections of chaotic and periodic 
oscillators, chaotic maps, etc) have recently been studied 
[lil HEl Hrl [I7L ITsj ] , and a transition to synchrony has 
been observed at coupling strengths that depend on the 
dynamics of the uncoupled systems. 

Our aim in this Letter is to generalize these previous 
works by studying large collections of heterogeneous gen- 
eral dynamical systems coupled by networks with com- 
plex topology. We find that in this general case there is a 
transition to coherence, and that the coupling strengths 
at which it occurs can be obtained from the uncoupled 
individual unit dynamics and the eigenvalues of the ad- 
jacency matrix. Thus, we achieve a separation of the 
problem into a part depending on the network only and 
a part depending on the individual unit dynamics only 
[l9| . Our model allows strong heterogeneity and different 
dynamics in the collection of dynamical systems, making 
it potentially appropriate to describe biological or other 
strongly heterogeneous systems. 
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We study networks of N coupled maps satisfying 

x% = f(x«\» j )+wW+ (1) 



m— 1 

Here j = 1,2, ... ,N labels the map and we are interested 
in large N, f(xn\ /ij) determines the uncoupled dynam- 
ics of each individual map with a parameter \ij ; for each 
map j the vector of parameters fij is chosen randomly 
and independently of j with a probability distribution 
p([i)', and the term iu„ is a random noise which is as- 
sumed to be statistically independent of j and n and to 
satisfy E[w^] = 0, E[w^w^}] = a 2 S nm Sji, where E[-} 
represents the expected value. In the coupling term, k 
is a global coupling strength and the scalar functions g 
and q are assumed to be smooth and bounded. The no- 
tation (• } represents the average over the distribution of 
the vector of parameters \ii and over the natural mea- 
sures of the attractors of the noisy uncoupled (k = 0) 
system. Alternatively, (q) is the /Lt-average of the infinite 
time average of q(x n ) over a typical orbit x n of the noisy 
uncoupled system. We remark that (q) is independent of 
time. We also note that x can be a vector for the situa- 
tion of multidimensional individual maps. For simplicity, 
in what follows we consider x to be a scalar. 

The N x N matrix Aj m determines the network of 
interactions: node m interacts with node j only if Aj m ^ 
0. We refer to those nodes m for which Aj m ^ as the 

neighbors of node j , and to dj = X)m=i A?™ as the degree 
of node j. 

We are interested in studying system Q for the case 
in which nodes have a large number of neighbors. In 
this case, if the initial conditions are chosen distributed 
according to the natural measure of the attractors of the 
uncoupled systems, then, because of the large number of 
terms in the sum in the coupling term in Eq. QJ, the 
fact that distributed according to the measure 

of the uncoupled attractors, and the lack of correlations 
between the parameter vectors and the network, we can 
approximate 

N N 

J2 A jm q(xir } ) « (q(x)) J2 A om, (2) 
m— 1 Tci— 1 

and equality of these sums applies in the limit of an in- 
finite number of neighbors. We refer to this situation as 
the incoherent state, and we will study in what follows 
its linear stability. Under the previously mentioned as- 
sumption of a large number of neighbors per node, the 
incoherent state is (approximately) a solution of the sys- 
tem Its linear stability can be studied using the 
same methods that were used in Ref. ^7| for the all-to- 
all case. In the following, we will adapt these techniques 
to the case of general connectivity. 



In order to study the linear stability of the incoherent 
state, we assume that Xn is in the incoherent state and 
introduce an infinitesimal perturbation Sxn ' ■ Lineariza- 
tion of Eq. produces 

N 

5x% = f'(x<P,H)5x® +kg(x^)Y / A l3 q'(x^)8x^. 

(3) 

In order to solve Eq. we consider (as in the variation 
of parameters method for differential equations) a pertur- 
bation e$ of the uncoupled system e^ +1 = f'(xn \ Mi) e n^ 
with e = 1. Defining = r)8x$ and 

assuming exponential growth, so that r„ = j^r) n , we 
obtain for large n [2(| 

In order to proceed further, we will again use the as- 
sumptions of large number of neighbors per node and 
statistical independence of the network and the vector 
of parameters. As we did in Eq. J5J, we approximate 
Eq. {1} by 

/ ™ n'(r {l) )f {l) a(r ( ' i) )nP- n - 1 \ N 
7 (m) = k I 1 ^n+l^n+lff^P )V \ ^ 

\p=0 € p+l I i=l 

(5) 

[If the sum over i in Eq. |0J is imagined as approximating 
N times the expected value of a product of two random 
variables, Eq. © approximates N times the product of 
the two expected values as suggested by our assumption 
of their independence.] 

Let Q{r]) be the average in Eq. jnj. With m = n — p, 
and letting n — > oo, 

Q(v) = n- 1 ( £ -^^gWi).9(z«- m )?T m ) ■ (6) 

\ n ^n— m+1 / 

This quantity depends only on the dynamics of the in- 
dividual uncoupled oscillators, and also results from the 
analysis of the globally coupled case 171. From Eq. 
we obtain vP' = k c \jQ{rf)u^\ where u^' and Xj denote, 
respectively, the eigenvectors of A and their correspond- 
ing eigenvalues. The onset of instability of the incoherent 
state corresponds to \rj\ = 1, or rj = e lu} . Thus, network 
mode j becomes unstable at a critical coupling strength 
satisfying 

1 = MiQCe*"), (7) 

where the critical frequency lo is found by 
Im[XjQ(e lu )] — 0, and Im denotes the imaginary 
part. We are interested in the solutions k c of Eq. Q 
with the smallest magnitude. Typically, they correspond 
to the mode associated to the eigenvalue of largest 
magnitude, which is usually real |2l| . In this case the 
critical frequency is found from Im[Q(e LU) )\ — and is 
independent of the network. 
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We remark that the definition and numerical compu- 
tation of Q(ri) is in general nontrivial, since for \t]\ = 1 
(necessary for investigation of the onset of coherence) 
the individual terms in the sum in Eq. © diverge with 
to for typical initial condition in chaotic maps. However, 
for large enough \r]\ these terms decay exponentially with 
increasing m and consequently the sum and average can 
be interchanged, so that 



Q(rj) = if 



E 

m=0 



-q'(x n+ i)g( 



(8) 



j ^n — m+1 

In Refs. 0, 0| it is argued that the averages in the 
summand of JSJ decay exponentially with m, so that © 
can be analytically continued to \rj\ — 1 as desired. 

We note that the function Q(rj) depends exclusively 
on the dynamics of the uncoupled oscillators and their 
parameter distribution function p(p), while the eigenval- 
ues Xj depend exclusively on the network. Therefore, 
an independent treatment of these problems allows de- 
termination of the critical coupling strengths for the full 
system given by Eq. Q). 

We now present two examples illustrating our theory. 
In order to quantify the coherence, we define an order 
parameter r by 



Em=l c M'?(>™ m) ) 



<?(*)> 



N -} 




(9) 



where <• >t denotes a time average and the in-degree dj 
is defined by d m = EjLi Arm- Note that the numera- 
tor can be written as J2jLi El=i ^-jm[<l( x n) ~ (l{ x ))h 
and, therefore, r measures the rms of the coupling term 
in Eq. Q [aside from the factor g(x^')]. Thus, the in- 
coherent state corresponds to r « 0. We will investigate 
what happens to r as the coupling strength k is increased 
past the critical values predicted by the theory. 

In our numerical experiments, we compute the order 
parameter r using Eq. © with x^ obtained from it- 
eration of Eq. We calculate the time average us- 
ing 1000 iterations after the initial transients have dis- 
appeared. As an example, we consider ^3] f° r the func- 
tions in Eq. f(xn,m) — 2x$ q(x) = cosx, and 
g(x) = sin(2a;)-|-sin(4a;). Throughout x is regarded as an 
angle-like variable and its value modulo 2ir to be taken 
where appropriate. In Ref. Q(e l ") was calculated for 
Gaussian noise with mean zero and standard deviation a 
to be 



Q(en = -- 



-a 2 /2 



-(cos(/j)) 



-5<t 2 /2 



(cos(3/j)) 



( 10 ) 

We numerically consider the following two exam- 
ples, listed below with the theoretical critical coupling 
strengths obtained using expression l|10|l in Eq. J7J). 

1. Identical noisy maps with a = 0.4, = <5(m)- 

For this example the solutions of Eq. are k c Xj w 




1.49 and fc c A,- 



-0.88. 



FIG. 1: The order parameter, r 2 , plotted on a linear scale 
as a function of the coupling strength k. The insets show the 
same data plotted on a logarithmic scale for r 2 . (a) Exam- 
ple 1 (identical noisy maps) with a scale-free network with 
N — 10 5 , exponent —3 and d min — 100. (b) Example 2 
(heterogeneous noiseless maps) with a scale-free network with 
N = 20000, exponent —2.5, and and lower cutoff d m i n — 50 
(boxes) and d m in = 200 (stars). The vertical lines indicate 
the theoretical values for the critical coupling strength. 

2. Heterogeneous noiseless maps with a — 0, p(fj) = 
2/ir if < n < tt/2, otherwise. In this example 
we obtain k c \j = 3tt/5 and k c Xj = — 3tt/2. 

For the network connectivity, we consider a scale free 
network with exponent 7, i.e., a network in which the 
degree distribution P(d) satisfies P(d) cx g?~ 7 for d > 
d m i n and otherwise. We impose a lower cutoff d > d m i n 
so that our assumption of a large number of neighbors per 
node is satisfied. We use 7 = 3, d m i n — 100 for example 
1 and 7 = 2.5, d m i n = 50 and d m i n — 200 for example 2. 
In order to construct such networks we use the Random 
Graph model of Chung et al. j2^ . 

For scale free networks as described above, the largest 
positive eigenvalue A+ is significantly larger than the 
magnitude of the most negative eigenvalue A_, and so 
the mode desynchronizing first is the one associated to 
the largest positive eigenvalue. Consequently, the value 
of Xj used to determine the critical coupling strengths k c 
for these examples will be A+. 

In Figs, ^a) and (b) we show the square of the order 
parameter as a function of k for examples 1 and 2 with 
N = 10 5 and N — 20000, respectively (the insets show 
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the same quantity on a logarithmic scale). For example 
1 [Fig. I^a)] the order parameter grows for values of the 
coupling strength close to the positive and negative crit- 
ical values predicted by the theory (vertical lines in the 
figures). For example 2 [Fig.QJb)] the transition on the 
positive side is quite sharp and occurs very close to the 
theoretical value k c \ ~ 1.88, while on the negative side 
the transition, although somewhat less well defined, also 
occurs close to the theoretical critical coupling strength 
k c2 ~ —4.71. On the negative side, the transition is not 
so sharp. However, we observe that as we increase d m i n 
from 50 (boxes) to 200 (stars), the transition becomes 
sharper. Generally, we find that the agreement with the 
theoretical results improves as the minimum degree d m in 
and N become larger |2^. The plots in Fig. were pro- 
duced by starting in the incoherent state and increasing 
the magnitude of the coupling strength k from zero. We 
have found that the transitions in our examples are not 
hysteretic, i.e., the same behavior is observed if we use 
the synchronized state as an initial condition and de- 
crease the magnitude of the coupling strength. We note 
that in the globally coupled case, it was found that the 
computation of Q{rf) fails to converge for an ensemble 
of identical noiseless logistic maps [lTj. It was argued 
that this results from structural instability of the map 
and singularities in its invariant density, which make a 
perturbation approach questionable. Since the definition 
and numerical determination of Q(i]) in our case and for 



the globally coupled case are identical, the lack of conver- 
gence observed for this example in the globally coupled 
case will also occur in the case of a network. However, 
we note that a small amount of either noise or parameter 
heterogeneity was shown in 01 to restore the validity of 
the results. 

In summary, we have studied the onset of synchro- 
nization in large networks of coupled maps (the case of 
coupled continuous time oscillators can also be treated by 
these methods and will be discussed elsewhere) . We have 
found that the critical coupling strength at which the 
transition to synchrony takes place depends separately 
on the dynamics of the individual uncoupled oscillators 
and on the largest eigenvalue of the adjacency matrix of 
the network. Thus, we have achieved a separation of the 
problem of the stability of the incoherent state for net- 
works of coupled dynamical systems into a part depend- 
ing on the dynamics of the uncoupled individual units 
and a part depending exclusively on the network |l9j| . 
Our theory directly generalizes the Kuramoto model of 
equal strength, all-to-all coupled phase oscillators to the 
case of oscillators with more realistic dynamics coupled 
in a potentially complex network. The results we obtain 
suggest that knowledge of network properties that favor 
larger maximum eigenvalues can be used to promote syn- 
chronism in large networks of coupled dynamical systems. 

This work was sponsored by ONR (Physics) and by 
NSF (DMS 0104087 and PHY 0456240). 



[1] M. E. J. Newman, SIAM Review 45, 167 (2003). 
[2] A. -L. Barabasi, and R. Albert, Rev. Mod. Phys. 74, 47 
(2002). 

[3] A. Pikovsky, M.G. Rosenblum, and J. Kurths, Syn- 
chronization: A universal concept in nonlinear sciences, 
(Cambridge University Press, Cambridge, 2001). 

[4] E. Mosekilde, Y. Maistrenko, and D. Postnov, Chaotic 
Synchronization: Applications to Living Systems (World 
Scientific, Singapore, 2002). 

[5] Y. Kuramoto, Chemical Oscillations, Waves, and Turbu- 
lence, (Springer- Verlag, Berlin, 1984). 

[6] S. H. Strogatz, Physica D 143, 1 (2000); E. Ott, Chaos in 
Dynamical Systems, Second edition, Sec. 6.5 (Cambridge 
University Press, New York, 2002). 

[7] J. G. Restrepo, E. Ott, and B. R. Hunt, Phys. Rev. E 
71, 036151 (2005). 

[8] J. G. Restrepo, E. Ott, and B. R. Hunt, Chaos 16, 015107 
(2006). 

[9] T. Ichinomiya, Phys. Rev. E 70, 026116 (2004). 
[10] D.-S. Lee, Phys. Rev. E 72, 026208 (2005). 
[11] T. Ichinomiya, Phys. Rev. E 72, 016109 (2005). 
[12] A. Jadbabaie, N. Motee, and M. Barahona. Proceedings 

of the American Control Conference (ACC 2004). 
[13] Y. Moreno and A. E. Pacheco, Europhys. Lett. 68, 603 
(2004). 

[14] A. Pikovsky, M.G. Rosenblum, and J. Kurths, Europhys. 

Lett. 34, 165 (1996). 
[15] H. Sakaguchi, Phys. Rev. E 61, 7212 (2000). 
[16] D. Topaj, W.-H. Kye, and A. Pikovsky, Phys. Rev. Lett. 

87, 074101 (2001). 



[17] S.-J. Baek and E. Ott, Phys. Rev. E 69, 066210 (2004). 

[18] E. Ott, P. So, E. Barreto, and T. Antonsen, Physica D 
173, 29, (2002). 

[19] In this respect our result is in the same spirit as that of 
the 'master stability function' of L.M. Pecora and T.L. 
Carroll [Phys. Rev. Lett. 80, 2109 (1998)]. The main dif- 
ference is that we consider the onset of synchronism from 
an incoherent state of many potentially different chaotic 
or periodic systems, while Pecora and Carroll consider 
the stability of a fully synchronized state of identical os- 
cillators. 

[20] The derivation follows closely the all-to-all case in 
Ref. |l7j|. which assumes A nm = 1. Further details will 
be published elsewhere. 

[21] For symmetric matrices or matrices with positive en- 
tries the eigenvalue with largest magnitude is real. For 
a large class of nonsymmetric matrices with mixed pos- 
itive/negative entries, we have found that if the average 
of the nonzero elements of the matrix is large enough, 
the eigenvalue with largest magnitude is also real and it 
is well separated from the eigenvalue with next largest 
magnitude. For details, see the Appendix of Ref. 0. 

[22] F. Chung, L. Lu and V. Vu, Proc. Natl. Acad. Sci., 100, 
6313 (2003). 

[23] A quantitative treatment and theory for the effect of d m i n 
is given in Q for the case of the Kuramoto model, and 
we expect similar behavior when the dynamics of / at 
different parameter values are sufficiently similar. 



